Electronic zero-point oscillations in the strong-interaction limit 
of density functional theory 



00 

o 
o 

(N 
o 

Q 

m 



X3 
O 

o 



> 
o 

00 

O 

> 

X 



Paola Gori-Giorgi,^'^ Giovanni Vignale,^ and Michael Seidl'' 
^ Laboratoire de Chimie Theorique, CNRS, Universite Pierre et Mane Curie, 4 Place Jussieu, 75252 Paris, France 

^ Af deling Theoretische Chemie, Vrije Universiteit, 
De Boelelaan 1083, 1081 HV Amsterdam, The Netherlands 
"^Department of Physics and Astronomy, University of Missouri, Columbia, Missouri 65211, USA 
"^Institute of Theoretical Physics, University of Regenshurg, 93040 Regenshurg, Cermany 

(Dated: December 3, 2008) 

The exchange-correlation energy in Kohn-Sham density functional theory can be expressed exactly 
in terms of the change in the expectation of the electron-electron repulsion operator when, in the 
many-electron hamiltonian, this same operator is multiplied by a real parameter A varying between 
(Kohn-Sham system) and 1 (physical system). In this process, usually called adiabatic connection, 
the one-electron density is kept fixed by a suitable local one-body potential. The strong-interaction 
limit of density functional theory, defined as the limit A — > oo, turns out to be, like the opposite non- 
interacting Kohn-Sham limit (A — > 0) mathematically simpler than the physical (A = 1) case, and can 
be used to build an approximate interpolation formula between A — > and A — > oo for the exchange- 
correlation energy. Here we extend the exact treatment of the \ ^ oo limit [Phys. Rev. A 75, 042511 
(2007)] to the next leading term, describing zero-point oscillations of strictly correlated electrons, 
with numerical examples for small spherical atoms. We also propose an improved approximate 
functional for the zero-point term and a revised interpolation formula for the exchange-correlation 
energy satisfying more exact constraints. 



I. INTRODUCTION 

Kohn-Sham (KS) density functional theory (DFT) 
[U, 0, II] is a very successful method for electronic struc- 
ture calculations thanks to its unique combination of 
low computational cost and remarkable accuracy. In 
the Kohn-Sham formalism, the ground-state energy of 
a many-electron system in a given external potential 
Vljxt = X]t=i "^extlri) is rewritten as a functional of the 
one-electron density p(r). 



where the non- interacting kinetic energy functional [p] 
is obtained by replacing Vcc with zero in Eq. ([2|). 



Ts[p\ = min(*lrl^') 



(6) 



and the Hartree functional U[p\ is the classical electro- 
static repulsion energy 



|r-r'| ■ 



(7) 



E[p]=F[p]+ / dVp(r)i.ext(r), 



(1) 



where 



F[p] - min(*ir + v;c|*), 



(2) 



with the operators (in Hartree atomic units e = m = h 
ao — I used throughout) 



T 



i=l 

1 ^ 1 - S,, 



2 \ri - r, 



(3) 
(4) 



In Eq. ^ the minimum search is carried over all anti- 
symmetric wavefunctions yielding a given density p [3]- 
The universal functional F[p] of Eq. ^ is further divided 
into 



F[p]^T,[p] + U[p] + EM 



(5) 



The only quantity that needs to be approximated is the 
functional for the exchange-correlation energy, Ey^c [p] , de- 
fined as the quantity needed to make Eq. ([5]) exact. The 
great success of KS DFT in solid state physics stems 
from the fact that even the simplest approximation for 
Exc[p], the local-density approximation (LDA), already 
gives remarkable results for basic properties of simple 
solids. A fundamental step forward to improve the solid- 
state physics results, and to spread the use of KS DFT 
into the quantum chemistry world, has been the advent of 
generalized gradient approximations (GGA), which are, 
to a large amount, due to the work of John P. Pcrdew 
and his coworkers [1, d, 0| • 

Despite its success in scientific areas now ranging from 
material science to biology, KS-DFT is far from being 
perfect, and a huge effort is put nowadays in trying to 
improve the approximations for i?xc[/o] (for recent re- 
views see, e.g., 'qJ). The focus of a large part of 
the scientific community working in this area has shifted 
from seeking explicit functionals of the density (like the 
GGA's), to implicit density functionals that construct 
the exchange-correlation energy from the KS orbitals. 
For example, predicted atomization energies of molecules 
have been improved by meta-GGA's (MGGA) 0, [H 
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which make use of the orbital kinetic energy density, by 
hybrid functionals (see, e.g., which mix a frac- 

tion of exact exchange with GGA exchange and corre- 
lation, and by range-separated hybrids, in which only 
exact long- or short-range exact exchange is used (see. 

The next step [l9!| towards higher accuracy could be 
fully non-local functionals, which use 100% of exact ex- 
change (for a recent review, see [13 )■ Despite several 
attempts and the increasing understanding of the crucial 
problems [2l| , the construction of a fully non-local corre- 
lation energy functional compatible with exact exchange 
is still an issue. A possible way to address this problem is 
to use the information contained in the strong-interaction 
limit of DFT . To explain this strategy, we have first 
to recall an exact formula ^23,] for E^c [p] , 

S,e[p] = / d\Wx[p]. (8) 
Jo 

The integrand Wx[p] is given by 

Wx[p] = {^x[p]\V,,\^x[p])~U[p], (9) 

where ^'a[/o], for a given value of A > 0, is the wave- 
function that minimizes {'^\T + XVcc\'^) and yields the 
density p. If p is u-representable for all A > 0, ^a[p] is 
the ground state of a fictitious iV-clcctron system with 
the Hamiltonian 

Hx[p]=f + XV,, + V,iM (10) 
where the A-dependent external potential, 

N 

Kxt[p]-E^oxt([p];r.), (11) 

ensures that Hx[p] have the same given (A — 1) ground- 
state density p{r) for all A. When A = 0, the Hamil- 
tonian of Eq. (flOl) becomes the KS Hamiltonian, and 
z;^^°([p];r) = vks{^), the familiar KS potential, while 
for A = 1 we recover the Hamiltonian of the physical 
system. 

We can use perturbation theory to obtain an expansion 
of Wa[p] in powers of A starting from A = 0, 

Wx[p] = E^[p] + 2 XE^'^^ip] + 0(A2), (12) 

where E^[p] is the exchange energy and £^^^^[p]is the 
second-order correlation energy in Gorling-Levy [24* per- 
turbation theory. This perturbation series expansion, 
however, seems to have a finite radius of convergence 
(Ac) which for many atoms and molecules is less than 1, 
Ac < 1 ^,25^. Moreover, evaluating terms of ever higher 
order becomes impracticably expensive. Nevertheless, 
the exact lowest-order terms E^\p\ and E^^'^lp] can be 
used for an alternative approach [22l |. called interaction- 
strength interpolation (ISI), to approximate the inte- 
grand in Eq. The basic idea of ISI is to combine 



the A ^ limit of Eq. (|12p with the information from 
the opposite strong-interaction limit, A — *■ oo, to con- 
struct an interpolation formula for Wx [p] . This way, the 
information on the physical system at A = 1 is extracted 
from an interpolation between A — *■ and X —^ oo. 

In the strong- interaction limit, A — )■ oo, we will show 
in the next sections that M^a[p] has the asymptotic ex- 
pansion 

Wx[p] = Woo[p] + + 0{X-n, (13) 

where p > |. The expansion (1131) was justified from 
physical arguments in Refs. [26l. l27l|. and a simple ap- 
proximation for the two functionals Woo[p] and ^^^^[p], 
the point-charge plus continuum (PC) model [11], has 
been used for the ISI, yielding atomization energies with 
errors within 4.3 kcal/mol [22]. In a recent paper [29| . 
the functional Woo [p] of Eq. (fT3|) has been constructed 
exactly. The main object of the present work is the ex- 
tension of the exact treatment of Ref. [l^ to the next 
term, VF^[p]. 

The paper is organized as follows. In the next Sec. HIl 
we briefly review the results of Ref. [l^, recalling that 
the strong-interaction limit of DFT reduces to a 3iV- 
dimension classical equilibrium problem whose minimum 
is degenerate over a three-dimensional subspace. In 
Sees, mil and IIVI we define local curvilinear coordinates 
based on the local normal modes around the degenerate 
minimum. These local curvilinear coordinates will be 
used, in Sec. |Vl to expand the Hamiltonian of Eq. (fTU]) 
for A — > oo, up to the order A^/*. The corresponding ex- 
pansion of ^'a[p] is carried out in Sec. IVIl and the exact 
expression for [p] is obtained in Sec. IVIH where we 
also report numerical results for small spherical atoms, 
and we propose an improved PC functional for W(^[p]. 
In Sec. IVIIII we revise the interpolation formula for the 
ISI functional in order to satisfy the exact expansion of 
Eq. ini) up to 0(A"i). The last Sec. US is devoted to 
conclusions and perspectives. More details of the deriva- 
tion of our expansion are given in Appendix [X] and a 
fully analytic example is reported in Appendix IB] 



II. STRICTLY CORRELATED ELECTRONS 
(SCE) 

In the A ^ oo limit it has been shown [2^, that, 
in order to keep the N electrons in the given density p, 
the external potential in Eq. (jlOp must compensate the 
infinitely strong interelectronic repulsion, thus becoming 
proportional to A, 

!Mdlil=„g^^([p],r), (14) 

A — ^cxD A 

with a smooth finite function i'sce([p]j r). (For brevity, 
the argument [p] will be often dropped in the following). 
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The leading term in Eq. (jlOp when A ^ oo is then a 
purely multiplicative potential-energy operator, 



III. THE SCE POTENTIAL-ENERGY 
MINIMUM 



i?A^oo[p] = A(T4c + Vsce) + 0(VA). (15) 



The square |5'> 



of the corresponding wavefunc- 



tion is a distribution that is zero everywhere except for 
electronic configurations for which V^o + ^sce has its 
global minimum. In order to guarantee a given smooth 
density p{r) in such a "classical" state,, this global min- 
imum must be degenerate over a three-dimensional sub- 
space of R^^ [2§| (otherwise, the density would be a sum 
of delta peaks centered in the equilibrium positions of the 
N electrons). We call this classical state with a smooth 
density "strictly correlated electrons" (SCE). The square 
of the SCE wavefunction I^'sceMP = | limA^cc- *a[p] P 
reads 



|vI/scE(ri,...,r^)P = -iy^ / ds^J(ri-/p(i)(s)) 



X(5(r2 - /p(2)(s))...(5(rAr - fp(jv)(s)). 



(16) 



where fi, .., are "co- motion functions" , with fi(r) ~ r, 
and P denotes a permutation of {1, ...N}. This means 
that the N points ri,...,rjv in 3D space found upon si- 
multaneous measurement of the N electronic positions in 
the SCE state always obey the — 1 relations 



= f,(ri) {i^2,...,N). 



(17) 



If the — 1 CO- motion functions fi(s) satisfy the differ- 
ential equation 



p(f,(r))dV.(r)-p(r)dV, 



(18) 



together with special transformation properties [231 (see 
also Ref. [S^l), the SCE wavefunction of Eq. (fTB]) yields 
the given density p(r). One has then to find the initial 
conditions for the integration of Eqs. psp that minimize 
the expectation of Vcc- The leading coefficient IVoo[p] in 
Eq. (|13p has a simple analytic expression in terms of the 
fi(s) [see Eq. (ITS)) ], and has been evaluated for spherical 
atoms with up to iV = 10 electrons [29| . 

In order to treat the next leading term, W4>[p] 
Eq. Uni), we have to consider the next terms in the 
A — > oo expansion of the Hamiltonian of Eq. (jlOp . i.e., 
the kinetic energy T and the next orders of V^tt- Physi- 
cally, we expect that [p] is determined by zero-point 
oscillations around the degenerate SCE minimum. In the 
following, we give a formal justification to this physical 
argument. 



Writing r = (ri, ...jTn) E R"^^ = SI, we consider the 
asymptotic potential-energy function (/? — > R), 



Epotir) 



hm^ 

A — >oo A 



N 



1 - 6,. 



1 

2 ^ |r, - r, 

Ko + VsCE- 



E 

i=l 



i'scE(ri) 



(19) 



As said, the SCE external potential wscE(r) has the very 
special property that the function -Epot(r) has a degen- 
erate minimum i?scE on the 3D subset 



/?o - {/(s) I s e R'} c /2, 



(20) 



where /(s) = (s, f2(s), fAr(s)), with the R-^ ^ co- 
motion functions fi(s). In other words, for all r g /?0i 
the function Spot (n) assumes the same constant value 



SscE = W^\p\ + U\p\ + ;^«scE(f^(s)) (21) 



which, in particular, is its global minimum within i7. 
For illustration, see the analytical example of Eq. (jB5[) 
in Appendix [Bl 

In the very limit A — > oo, when ii\\p\ Ai?pot(zi) + 
0(-\/A), the square of the wave function |^'a[p]P becomes 
the distribution I^I/sceWP of Eq. p^ . which is strictly 
zero almost everywhere in Q except for the 3D subset J?o 
where i?pot(zi) is minimum [l^. 



«'sce([p],i:) = V7:ef?\/2o. 



(22) 



For large, but finite A 1, the electrons are expected 
to perform small zero-point oscillations about the SCE 
configurations r G J?o , within a narrow 37V-D "envelope" 
(with a small width e > 0) of the 3D subset Qq C i7, 



12, = {r G n\d(r_,QQ) < e}. 

Here, for a given r E il, the quantity 

d{r, Qo) := min |r- /(s)| 
seR3 — 



(23) 



(24) 



is the minimum 3iV-D distance between r and any /(s) G 
Hq. Notice that J7o C /2e C i7 and Hq — lim^^o i?,. 
For rGil^, i?pot(z:) may be expanded about r(s) G ilo, 



^ 3N 

x(r^-/p(s))(r,-/,(s)) + ... (25) 

Since Epot{r) is minimum at r = /(s), there are no first- 
order terms. [The dots represent the terms of third and 
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higher orders.] For any given s G R^, the Hessian matrix 
M^i^(s) in the second-order term has 3N non- negative 
eigenvalues a;^(s)^ which can be labeled such that 



a;^(s)2 = (a* = 1,2,3), 
u;^(s)2>0 (Ai = 4,...,3iV). 



(26) 



The corresponding 3A^-D normalized eigenvectors e''(s), 
with components e^{s) {a = 1,...,3N), are pairwise or- 
thogonal. 



3N 



e^(s).e^(s)^^e^^(s)<(s) = <5^.. (27) 

CT = 1 

The first three eigenvectors, with zero eigenvalues, lie 
in the space "tangential" to i7o, the remaining 3A^ — 3 
eigenvectors are "orthogonal" to i?o, 

e^(s) • = (m = 4, 37V, a = 1, 2, 3), (28) 

OSq. 

where a ~ 1,2,3 denotes the three cartesian components 
{x,y,z) of s. 



IV. LOCAL NORMAL MODES 

For sufhciently small e > 0, we use these eigenvectors 
to introduce a set of 3A'^ curvilinear coordinates in fi^. 
A given point r = (rn, ri2, ri3, tati, r-Ara, tats) € ^e, 
is written in terms of these local curvilinear coordinates 
as follows. The first three curvilinear coordinates are the 
cartesian coordinates si, S2, S3 of the minimizing vector s 
in Eq. p4|) . fixed by the condition that the 37V-D vector 
r — /(s) in n is orthogonal to i7o in the point /(s), 



(r-m) 



dm 

dSa 



(a = 1,2,3). (29) 



The remaining 3A'^ — 3 coordinates are the projec- 
tions qi,...,qzN of r — /(s) onto the local eigenvectors 
e4(s),...,e3^(s), 



.-/(s) = ^q^e^(s). 

P=4 



(30) 



The first three eigenvectors e^'^''^(s) are not needed, since 
they are tangential to Qq at the point /(s) and therefore 
orthogonal to r — /(s). Inverting Eq. J^Ol) yields 

q>.=d'-{r-l{s)) (Ai = 4,...,3iV). (31) 

For these new curvilinear coordinates, we also write 

(si,S2,S3,94,---,g3Ar) = (s,?). (32) 

Notice that r has 3A^ components, while q has only 3A^— 3 
ones. In this notation, Eq. (|5D|) reads 

3JV 

7v = /.(s) + ^ei;(s)q^ (i. = l,...,3iV). (33) 

p=4 



This is the transformation formula between the cartesian 
coordinates r and the "local normal modes" (s, q) in the 
3iV-D configuration space Q. 

In terms of the g^, the second-order contribution in the 
Taylor expansion ((25|) becomes diagonal, 

3Af 



M=4 



^ ZN 



(34) 



Here, the third-order term is derived from the corre- 
sponding term in Eq. ()25p (in the present notation). 



^ 37V 



53i?pot(r) 



3! 



^1 dr^drridr^ 



r=f{s) 



x(rc - /c(s))(r, - /,(s))(rc - /c(s)). (35) 
Using here Eq. ([33|) for Vi, — f^{s), we find 

3N 



e^^(s)e:;(s)e^(s). 



(36) 

Substituting Eq. (l33|) for r in the wave function ^'a(z:) 
that represents the state 4'a[p] yields the transformed 
wave function ^'a(s, q). While the original wave function 
obeys 



dVi... / dVw|*A(r) 



dr\-^x{r)\' 



the transformed one is normalized according to 



d^s / dqj{s,q) |*a(s, g)| = 1 



(37) 



(38) 



where J(s, q) is the Jacobian associated with the coor- 
dinate transformation ([33|) . see Eq. (IA12|) in Appendix 

El 

For sufficiently large A 3> 1, the wave function 5'a(s, q) 
strongly suppresses all configurations r € O except for 
the ones inside the narrow envelope i7e of the 3D subset 
Hq. This means that vI'a(s, q) is essentially different from 
zero only for (q^ + ... + q^j^Y^^ < e, where e decreases 
with growing A 3> 1 and goes to zero in the limit A ^ oo. 

More precisely, since the quadratic term in Eq. ([M]) 
is multiplied by A in the Hamiltonian (fTO)) . the scale of 
the quantum fluctuation is e ^ A^^/* = a for A — > oo. 
Therefore, it will be useful to switch for a given value of 
A 3> 1 from the present curvilinear coordinates (s, q) to 
scaled coordinates (s, u) where 

u = X^/^q ^ q = au (a = A^^/^). (39) 

This second transformation yields the wave function 

W„(s,u) = $A(s,aw). (40) 
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According to Eq. p8|) . we now have 

d^s J duKa{s,u)\^a{s,u)\'^ ^ I, (41) 

with the scaled Jacobian 

Kais,u) = a^'^-^J{s,au). (42) 
Later on, we shall make use of the expansion 

37V 

J(s,g) = J(s,0) + J2 ^^'^8)9^ + 0{ql), (43) 
whose derivation is reported in Appendix [Xl 

V. EXPANSION OF THE HAMILTONIAN 

To obtain an expansion for large A 3> 1 (or, equiva- 
lently, for small a = A^^/^ <C 1), we must express the 
Hamiltonian H\ [p] of Eq. (|10p in terms of the scaled co- 
ordinates (s,u). To this end, we split Hx[p] into three 
pieces, 

Hx [p]=f + XEpot (r) + (Kxt - A VscE ) , (44) 
and treat these separately now. 

A. Kinetic energy (first term) 

For the kinetic-energy operator T, the 3iV-D Laplacian 
is obtained in Appendix A in terms of the curvilinear 
coordinates from the general transformation rule 



3 3^ 



■AN 



1 d 



d 



(45) 

(To simplify the notation, we write s^i = for p = 1,2,3 
in this subsection.) Here, the matrix G'^" is the inverse 
of the metric tensor G^j/, defined by 

«--eI^|^-|^-|^- (46) 

^ dg^ dq^ dq^ dq„ 

and G is its determinant, G = det(G^j,). Switching in 
a second step from the q^ to the scaled coordinates 



yields the expansion (see Appendix [X)) 

f = VA[f + af^'^ + a^f^^^ + 0(a3) 
The operators T*^"-* are independent of A (or a 



f{Q) ^ _i.V — 
fj.=4 



f (1) 



1 Pi 



(47) 
A-V4), 

(48) 
(49) 



where ^^(s) is reported in Appendix lAl Notice that the 
a term is constant, since a 



B. SCE potential energy (second term) 



For the second term in Eq. (|44p . we use the Taylor 
expansion (|34p . with q^ = au,,. to find 



2 3N 



SCE ■ 



E'^^( 

M=4 



8)2^2 



3N 



4 3Ar 



+ ir E 4tlr(sK^.W + 0(a5) . (50) 



C. The remaining external potential (third term) 



For the last term in Eq. ((44|) . we make an ansatz that 
will later on turn out to be consistent. 



V;^.t-AVscE = A/A^a"\^(")(r:). 



(51) 



ri=0 



Using Eq. pop for r and q^ = au^, we may expand 



3N 



V(")(r) = V"^"'(/(s)-f^a^e'^(s)i 

3N 

= l^("H/(s)) + a5]V^i")(/(s))^e;^(s)t 



3iV 



a=l fi=4 
2 3iV 3iV 

E ^i"H/(s)) E e(^(s)e^(s)z.,^.. 

(T.T—l fl,V — 4: 

+0{a^). 



(52) 

Here, the coefficients vi"\ vir\ etc. denote the partial 
derivatives of ^("^(r) at r = /(s), 

l/(;)(/(s)) = — -A=i etc. (53) 

- dr^drr £=/(s) 

Now, Eq. ((5T|) yields the expansion 
K-^t-AFscE - VA[]/("'+aF(i)+a2i/(2)+o(a3)] , (54) 

with a-independent (multiplicative) operators 
Vio) ^ ^(o)(/(s)), 

3N 3N 



yd 



CT = 1 

3N 



M=4 
3iV 



= i^(2)(/(s)) + ^yW(/(s))^e;^(s)u^-f 



CT = 1 



M=4 



^ 3JV 3A' 

- J2 V^r\m) E e'^i^Ki^W^^- (55) 



/j,;/=4 
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D. Full Hamiltonian 

Eventually, combining Eqs. (|47| . (|50)) . and (|54)) . we 
obtain the expansion (recall that a — A^^^"') 

Hx[p] — Ai?scE 

+ x/A[7?(")+aH(i)+a27j(2)+o(c,3)l (56) 

with a-independent operators H*-"'. The first two terms 
read 



3N 



/i=4 



= 



9 



4 i:^»<"a^ 

/I— 4 ^ 



3iV 



3W 



^yj")(/(s))5]e(^(s)^ 



^=4 



^ 3N 



3! 



(58) 



VI. EXPANSION OF THE GROUND STATE 

Due to Eq. ([TO]) , the lowest eigenvalue Ex[p] of -ffA[p] 
(i. e., its ground-state energy) has the expansion 

Ex[p] = A£;scE + 

+ \/A + aE<-^^ + a^E^^^ + 0{aH . (59) 

We define E'^ [p] = E^°^> + aE'^^^ + a^E^^^ + 0{a^) as the 
lowest eigenvalue of the operator 

Since E'sce is a constant, -ffA[p] and H'^[p], with a = 
A-i/4, have the same ground state 



*a(s,u) 



^(0) +q2^(2) ^Q(^3) 



(61) 



For the a-dependent normalization constant, 

TV; = y (fs J duK^{s,u) ¥°'^{s,u) + 0{a) ^ (62) 
we obtain 



AC. = 



[1 + 0{a)] 



when '9^^^ is normalized according to 

I du J(s,0) ■^'^°Hs,u) ^ = 1. 



Collecting terms of equal orders 0(a") in the eigen- 
value equation H'^[p]'^a — yields a hierarchy 
of equations. The leading one is 7f(o)*(o) = E^°H^°\ 
where i?'"-' is given by Eq. ((57|) . For a given fixed 
s £ R^, the Hamiltonian H^'^'> describes an uncoupled set 
of 37V — 3 harmonic oscillators in ID. To be more precise, 
these oscillators are coupled via the dynamical variable 
s, but the dynamics of s is much slower, only appear- 
ing at orders 0{X°). Consequently, the leading term in 
the wave function factorizes into a product of Gaussians 
^Uu) = (^)i/^e-"'/2, with JZ,du\<t>Uu)\^ = 1, 



3N 



f(°)(s,u)=c(°)(s)n<fc.,(s)(^ 

M=4 



(65) 



Since T^^"-* (/(s)) is a pure multiplicative operator, the 
resulting eigenvalue of i?'^*'^ is 



3N 



= y(o)(/(s)) + ^ 



(66) 



M=4 



Due to Eq. (|59p . this expression cannot depend on the 
variable s, implying a condition on the n = coefficient 
in our ansatz (|5ip . 



3N 



In particular, we have 



t^^(s) 



const VseR^. (67) 



i3 



N 



3N 



w^(s) 



^=4 



(68) 



The role of the external potential at the order -s/A in 
Eq. pO|) is thus to keep the degeneracy of the SCE mini- 
mum (found at the order A) through the order v'A. This 
is necessary in order to keep the given smooth density 
p(r): if one of the SCE configurations (i.e., a given par- 
ticular So) had a lower energy than the others, the SCE 
wavefunction would collapse in that particular Sg, and 
the density would become a sum of delta peaks centered 
in f,(so) (with i = 1,...,N). 

In order to determine the prefactor C'"^ (s) of the 
wave function (|65p we observe that in the wave func- 
tion ^a(s,(7), the coordinate s € has the probability 
distribution 



(63) 



(64) pa(s)= / dqj{s,q) 



Pa(s) = J J(s,2)|*A(s,g)P 

= / dqJ{s,q)\^a.{s,\'/\)\^ (69) 



where a = \ Using Eqs. ((6T|) and ((63)) . we find 

|*(")(s,Ai/4^)|2 



^3N-3 



l + 0{a) . (70) 



7 



In the hmit A — s- cxd when p\{s) must become rigorously 
proportional to the electron density p{s), 



hm px{s) = — 



(71) 



the terms 0{a) in Eq. ([70| can be dropped and Eq. ((65|) 
yields 

^ = hm |C(«)(s)p / dqj{s,q)x 

iV A-+00 J — — 

3N 



(72) 



Since ^ujiu) is a normalized Gaussian, the /i-th factor of 
the product in Eq. ((72| approaches the (5- function (^(g^) 
as A ^ oo. Therefore, the right-hand side of Eq. ([7^ 
equals |C(°^(s)p J(s, 0), implying the result 



|d°)(s)|^ 



1 P(s) 
TV J(s,0)' 



(73) 



The next order in the perturbative treatment of the 
ground-state energy of Eq. ((60)) leads to 



£;(!) ^ = T/(i)(/(s)). (74) 

The same argument used for Eq. (|67p yields 

F(^'(/(s)) = const., (75) 

independent on s. The important point here is that the 
terms coming from T and 14o in the Hamiltonian ij'^^ 
of Eq. have zero expectation on the ground-state of 
the harmonic oscillator, so that there is no contribution 
to this order to the large-A expansion of Wa[p]. As we 
shall see in the next Sec. IVIIl the order \/Aa = A^/** in 
E\[p] of Eq. ([59]) corresponds to the order X^^^^ in the 
large-A expansion of VFa[p]. 

Notice that, in our treatment of the strong-interaction 
limit of DFT, we did not consider the effect on the energy 
of the spin state or, more generally, of the statistics. This 
is because the electrons are always localized in different 
regions of space well separated from each other. The 
effect on the energy of the spin state or of statistics in 
the A — *■ oo limit can be estimated as being of the order 
©(e""^^^"), which is the order of magnitude of the overlap 
between two different gaussians of Eq. . 



VII. THE COEFFICIENT W^p] 

From the expansion of E\ [p] of the previous Sec. lVIi wc 
can easily compute W\ [p] using the Hellmann-Feynmann 
theorem: 



Wx[p] 



dEx[p] f , .dv. 



cxt(r)^3^_ 



dX 



(76) 



From Sec. IVIl we obtain, in the A — > cx) limit. 



Ex[p] - / p{r)v^^dr)d'r = A(vI'sce| V;e|*scE) + 



0(A") 



(77) 



By differentiating both sides with respect to A, from 
Eq. ([75)1 we obtain the expansion for W\[p] of Eq. 
with 



N N 



i=l jyi 



f,(s)l 



U[p], (78) 



in agreement with the results of Ref. [29] , and the exact 
expression for the next leading term. 



wL[p]=Ud^s(ff:^^^'^ 



(79) 



This result generalizes (and proves) Eq. (35) of Ref. [2y| 
for spherical two-electron densities. As shown by 
Eq. (ini), there is no A"^/"* term in Wx^oo[p\- There is 
also no term oc A^^, which would imply a term oc log(A) 
in E\[p] and thus in the kinetic energy {'^ \\T\^ \) . Such 
a term would violate the known high-density scaling of 
(^Alrl^-A) [aH (see also the erratum of Ref. (H). 
As an example of application, we have computed 
for the same set of spherical (or sphericalized) 
atomic densities used in Ref. [29|] to evaluate Woo [p] ■ For 
each point (fi(s), ...,f7v(s)) on the degenerate SCE min- 
imum constructed in Ref. [295, we have evaluated the 
hessian matrix, the corresponding eigenvalues (s) , and 
thus W'o^ [p] of Eq. ([75)1 . In Table|T]we compare our results 
with the approximate PC functional ^8| , 



'PC 



|Vp(r)| = 



p(r)7/6 



(80) 



where C = 1.535, D = -0.02558. 

As explained in Ref. [29j . the SCE minimum for spher- 
ical densities is constructed from a set of radial co-motion 
functions and the angular minimization is done numer- 
ically. When one of the electrons is close to the nu- 
cleus, the numerical minimization displays instabilities 
in the smallest (but non-zero) eigenvalues of the hessian. 
However, as shown by Eq. ([7^ . such configurations are 
weighted by the density (in the spherically symmetric 
case by 47r s^p(s)) so that the error they introduce is rel- 
atively small. This error, however, increases with the 
number of electrons. The number of digits in our results 
of Table [J is determined by this numerical error. Notice, 
however, that Table |T] shows that the discrepancy of the 
PC model with respect to our results is much larger than 
our estimated numerical errors on the SCE values. 

While the PC model for the coefficient M^ooM makes 
errors of the order of 60 mH [29|] , we see from Table U 
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SCE (H) PC (H) error (mH) 



haviors, 



He 0.62084 


0.729 


108 


Li 1.38 


1.622 


240 


Be 2.59 


2.928 


334 


B 4.2 


4.702 


502 


C 6.3 


7.089 


840 


Ne 22 


24.423 


2423 



TABLE L Comparison of the values H^4;[p] in Hartree atomic 
units obtained with the SCE construction, and with the PC 
model The absolute errors of the PC model are also 

reported. 



SCE (H) revPC (H) error (mH) 



Li 1.38 


1.4066 


26 


Be 2.59 


2.579 


11 


B 4.2 


4.207 


7 


C 6.3 


6.43 


130 


Ne 22 


22.96 


960 



TABLE II: Comparison of the values W4, [p] in Hartree atomic 
units obtained with the SCE construction, and with the re- 
vised PC model of Sec. lVHl The absolute errors of the revised 
PC model are also reported. 



that the functional W4j[p] is much more seriously over- 
estimated. We can reduce these errors by recalling that 
in the PC model for W^[p\ the coefficient D of Eq. jHO]) 
was fixed by the condition that the PC value for the 
He atom be equal to the one obtained from the MGGA 
functional of Ref. 10]. Now that we have exact values, 
it seems natural to change D in order to make the PC 
model equal to the SCE result for the He atom. This 
gives D = —0.028957. The values for the other atoms 
obtained with the revised PC model are reported in Ta- 
ble [Hi we see that the error is now substantially reduced. 



VIII. REVISED ISI 

In Refs. [l^, [2^ an expression for that inter- 

polates between the two limits of Eqs. (fT2| and (fT3)) 
has been proposed and tested using the PC approx- 
imation for the functional T4^oo[/o] and The 
interaction-strenght interpolation (ISI) formula for W\ [p] 
of Refs. [12, H^, however, contains a spurious term c>c 
in its A — > cxD expansion [2^ . which, as explained after 
Eq. ([75)l . has the wrong scaling behavior in the high- 
density limit. Here we propose a revised ISI functional 
which does not have this problem. 

Instead of modeling W\[p],we use the same ISI in- 
terpolation formula of Ref. [22| directly for the integral 



EM 



dX'Wx'ip], 



(81) 



AJSI 



[p] = a[p] A - 



b[p]\ 



y^l + c[p]X + d[p] ' 



(82) 



The four functionals a[p], b[p], c[p] and d[p] are deter- 
mined by imposing the A — > expansion of Eq. (|12p and 
the A — !■ 00 expansion of Eq. (|13p . and they are thus 
determined by the two weak-interaction limit functionals 
Ex[p] and E^^'^[p] and the two strong- interaction limit 
functionals VFoob] and W!^[p], 

a[p] 
b[p] 



Woo[p] 

8E^^'[p]W:^ 



d[p] 



iE,[p]~w^[p]r 

leE^^^lpfw^p]^ 
{EM-Wo,[p]r 

HE^^'[p]W:^[p]' 



-1 



iE4p]~Wo,[p]r- 



(83) 
(84) 

(85) 

(86) 



The final formula for the revised ISI functional is ob- 
tained by putting A = 1 in Eq. ((82)) . 



e: 



rcvISI 



b[p] 



y/1 + C[p] + d[p] 



(87) 



satisfying the exact A — > and A — s- cx) asymptotic be- 



For the correlation energy of the neutral atoms consid- 
ered here, this revised ISI gives essentially the same re- 
sults of the original ISI of Ref. 



IX. CONCLUSIONS AND PERSPECTIVES 

We have presented an exact treatment of the strong- 
interaction limit of density functional theory up to the 
second leading term, describing zero-point oscillations of 
strictly correlated electrons. We have evaluated numeri- 
cally this zero-point contribution for small atoms, and we 
have used our results to improve a previous approximate 
functional for this term. A new interpolation formula 
for the exchange-correlation energy, satisfying more ex- 
act constraints, has been proposed, and will be tested 
elsewhere. 

Besides the possibility of constructing an interpolation 
formula for E^c [p] , the two functionals Woo [p] of Ref. [2§| 
and [p] evaluated in this work, are of valuable interest 
for the development of Kohn-Sham DPT. They are an ex- 
ample of exact implicit density functionals for systems in 
which the electron-electron repulsion largely dominates 
over the kinetic energy. They can be used to test proper- 
ties of the exact exchange-correlation functional like the 
Lieb-Oxford bound [s^, 'ss'l , and to test how approximate 
functionals perform in this limit [34, 3^ . 

Several issues still need to be addressed and will be 
the object of future work. The main problem of the ISI 
functional is the lack of size consistency. In order to be 
size-consistent, the interpolation of Eq. ([82]) should be 
done locally, using energy densities all defined in the same 
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gauge. A first step in our future work will be the analysis 
of exact energy densities for the functionals Woo [p] and 
(see also Ref. [IBl), and the construction of corre- 
sponding approximations. Another important problem is 
the development of a reliable algorithm to solve the SCE 
problem for a given non-spherical density. Other promis- 
ing research lines are the study of the next leading term, 
which is of purely kinetic origin, and the construction of 
approximations to describe the effect of the spin state on 
the energy. 
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APPENDIX A: TRANSFORMATION OF THE 
LAPLACIAN 

In order to write down the components of the metric 
tensor dk of our local curvilinear coordinate transfor- 
mation, we define the indices as follows: a, (3, 7, .. denote 
the cartesian components 1,2,3 = x, y, z of s, the indices 

v,a,T, ... denote the normal-mode components q^, and 
the latin indices k, ... denote general components, either 
a, .. or /i, .... We then have to distinguish three blocks in 
the metric tensor Gik'. af), p,u, and a/j,. 



3N 



M=4 



ds^dsp 



e'^(s) + 



3Af 

E 

3N 



M=4 



(Al) 

(A2) 
(A3) 



where in Eq. (|Aip we have defined the 3x3 metric tensor 
(7q,/3(s) which only concerns the coordinates si, 52,53, 



5a/3(s) = 



dm dm 

dsa dsp 



(A4) 



When A — > cx), our wavefunction is zero everywhere ex- 
cept very close to J?Oj i-e., for very small cx A^^^^. 
Introducing the scaled coordinates = X^^'^q^, we see 
that the metric tensor Gik has the A-dependence 

Gik = Gf^ + — ^ WM^ifc + TT/J E ''^l-'^'^^ik ^ 



M=4 



(A5) 



where and Z^" are tensors of elements 



and 



_ 9em .. . 



(A6) 

(A7) 
(A8) 

(A9) 

(AlO) 
(All) 



and G^^^ has elements G^'j, = gap, G^°j = and all 
the off-diagonal components equal to zero. In order to 
compute the large-A expansion of Eq. ([i5|) . we have to 
expand the determinant G, and the components G^^ of 
the inverse metric tensor. Using standard formulas, we 
obtain 

(1 \ 
Ai=4 Q/3 J 

(A12) 

where g is the determinant of gap, and g"^ are the com- 
ponents of its inverse. The tensor of components 
G*'^ has the large-A expansion, up to orders \~'-^/^, 



G-i = G('')"'--ij^^/^G(°) 'a^G^") \ (A13) 



3N 



Ai/4 



p=4 



Inserting these expansions into Eq. ([15|) we obtain 
Eqs. dlH]) and (gH]) with 

X'^{s) = ^J29"''i^)KM (A14) 

a/3 

Finally, the Jacobian of our change of coordinates is 
simply equal to VG of Eq. (|A12p . 



APPENDIX B: ANALYTIC EXAMPLE 

As an illustration, we consider a system of two elec- 
trons in ID space (i.e., on the cc-axis) with a given 
ground-state density p{x), 



dx p(x) — 2. 



(Bl) 



In this case, Eq. p?)) reads X2 ~ f2{xi), with the single 
co-motion function /2(s) = f{s) which, according to Ref. 
[29! . obeys the differential equation p(/(s))/'(s) = p{s). 
For the Lorentzian density, /(s) is found analytically. 



p{x) 



2 1 

n I + x^ 



fis) 



(B2) 
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In this case, the SCE external potential, fixed by the con- 
ditions 3|wscE(a;) = sgn(a;)|x-/(x)|-2 and vsce{x) 
for a; — > ±00, is given by 



vsce{x) — arctan(a;) — 



1 



TT 

2' 



(B3) 



In terms of /(s) = (s,/(s)), Eq. (^0)1 now yields a ID 
subset oi Q = R^, 



= {/(s)|s e R} C 



(B4) 



In the example (IB2p . J7o is given by the two branches of 
the hyperbola X2 = f{xi) = — ^ in the a;iX2-plane il. 
In the following, we focus on the branch Hq with xi > 
and X2 < 0, = {/(s)|s G R+}- 

The asymptotic potential-energy function, cf. Eq. 



+ t'scE(a;i) -f wscE(a;2), (B5) 



assumes its highly degenerate minimum for all x G i?o- 
Consequently, the first partial derivatives. 



dxi 
dx2 



(xi-.X2)2 (l + a;2)2' 
1 



x\ 



2\2 ' 



(B6) 



are vanishing for x — f{s) when the Hessian matrix of 
Epot[x) becomes 



It has the two eigenvalues 

= 0, U2{sf 



1 



2s 



(B7) 



(l + s'^)>0. (B8) 



(l + s2)3 

The corresponding normalized eigenvectors are 
1 



s 

A\ 1 



(B!)) 

While e^(s) is tangential, e^(s) is orthogonal to J7q^ at 
f{s)£ (2q and generally given by 



(BIO) 



For a point a; = (a;i,X2) G ^e, close to the curvi- 
linear coordinates (s,g) are defined by Eq. ([SO]) . 



/(s) + e(s)g 



(Bll) 



where s is fixed by the condition that the vector e(s) in 
the a;iX2-plane is orthogonal to at f{s) S . 

In terms of the partial derivatives of Eq. (jBlip . the 
metric tensor is given by the (2 x 2)-matrix 



/ dx dx dx dx 

ds ds ds dq 
dx dx dx dx 



(B12) 



V dq ds dq dq 
Using Eqs. (IbTO]) and (|BTT|) . we obtain 



dx 
ds 



^ = /'(.s) + ge'Gs) 



dx 
dq 



1 + 91 



1 + /'(s)2]3/2 



and thus 



1 



(B13) 



(B14) 



where g{s, q) = J{s, g)^, with the Jacobian 



J(s,g) = 



1 + 97 



v/T+iW- (B15) 



In the particular case of the density (|B2[) . we have 

2s 



J(s,g) 



1 



(B16) 



and the coefficients of Eq. (|T3|) are given by 

p(s) 1 



/•OO 

W^ooW + t/W = 2 / ds 
Jo 



2 s-/(s) 

1 f°° ds2s 1 



0.31831. 



dsP^UJ2is) 



TT 



(B17) 







2 

2 f"" ds 
^Jo JlTs^ 
0.633902. 



2s- 



l + s^ 



(BIS) 
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